clear all;
load A;
load b;
load c;
eqn.A_=A;
eqn.B=b;
eqn.C=c;
load N;
eqn.N_=N;
opts.blyap.maxiter=10;
% % Set tolerances and iteration limit for the Lyapunov-plus-positive solver
% opts.blyap.res_tol  =  1e-6;
% opts.blyap.rel_diff_tol  =  1e-6;
% opts.blyap.maxiter = 6;
% 
% % Set options for the ADI-based inner Lyapunov solver
% % and related shift computation
% opts.shifts.num_Ritz = 50;
% opts.shifts.num_hRitz = 25;
% opts.shifts.method = 'projection';
% opts.shifts.num_desired = 6;
% 
% opts.adi.maxiter = 8;
% opts.adi.res_tol = 1e-12;
% opts.adi.rel_diff_tol = 0;
% opts.adi.info = 2;
% opts.adi.norm = 'fro';
% opts.norm = 'fro';
% 
% % Set mess_res2_norm options
% opts.resopts.res.maxiter = 8;
% opts.resopts.res.tol = 1e-6;
% opts.resopts.res.orth = false;
% 
% opts.fopts.LDL_T = false;
% opts.fopts.norm = 'fro';
% 
% opts.srm.tol = 1e-10;
% opts.srm.info = 1;
[oper, opts] = operatormanager(opts, 'default');
[out,~,~,~]=mess_lyapunov_bilinear(eqn,opts,oper);
%[out,~,~,~]=mess_lyapunov_bilinear(eqn);